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ABSTRACT 

Methods  of  accounting  for  load  variation  in  vibration  signals  from  helicopter 
transmission  systems  are  presented.  These  methods  are  based  on  autoregreassive 
moving-average  (ARMA)  models,  and  several  ARMA  parameter  estimation  schemes 
are  presented.  Simulations  of  load  variation  are  carried  out,  and  a  prediction  error 
filter,  based  on  the  ARMA  models,  is  used  to  generate  a  residual  signal.  Fault  indices 
extracted  from  the  residential  signal  are  used  to  indicate  the  presence  or  absence  of  a 
fault.  The  results  of  the  simulations  suggest  that  this  method  of  fault  detection  is  able 
to  detect  both  general  and  local  fault  conditions. 


RELEASE  LIMITATION 

Approved  for  public  release 


Department  of  Defence 
- 4. - 

DEFENCE  SCIENCE  AND  TECHNOLOGY  ORGANISATION 


Published  by 


DSTO  Aeronautical  and  Maritime  Research  Laboratory 
PO  Box  4331 
Melbourne  Victoria  3001 

Telephone:  (03)  9626  7000 
Fax:  (03)9626  7999 
©  Commonwealth  of  Australia  1997 
AR-010-104 
February  1997 


APPROVED  FOR  PUBLIC  RELEASE 


Signal  Processing  Methods  for  Gearbox  Fault 
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Executive  Summary 

Current  vibration-based  techniques  used  for  gearbox  fault  detection  rely  on  analysis 
of  signals  captured  under  fixed  load  and  fixed  speed  conditions.  While  helicopter 
gearboxes  operate  nominally  at  fixed  speed,  they  do  not  operate  under  fixed  load.  If 
such  fault  detection  techniques  are  to  be  implemented  as  part  of  self-contained  on¬ 
board  systems,  then  adaption  to,  or  accounting  for,  load  variation  is  a  necessary 
requirement.  This  report  details  the  signal  processing  methods  which  are  able  to 
account  for  load  variation. 

Several  methods  of  linear  modelling  are  introduced.  These  techniques  enable  the 
vibration  signals  recorded  from  helicopter  transmission  systems  to  be  modelled.  The 
aim  is  to  identify  a  suitable  model  of  the  gearbox  vibration  signals  under  normal 
operating  conditions,  which  is  able  to  be  adjusted  to  meet  variations  of  the  load 
imposed  on  the  gearbox.  Once  a  suitable  model  is  established,  the  vibration  signal 
recorded  from  the  gearbox  is  compared  with  the  signal  produced  by  the  model. 
Various  signal  characteristics  (referred  to  as  fault  indices)  are  derived  from  this 
comparison  as  a  means  of  deciding  if  a  fault  is  present  or  not. 

The  results  of  several  simulations  show  that  this  type  of  linear  modelling  provides  a 
useful  means  of  both  accounting  for  load  variation,  and  detecting  the  presence  of 
faults. 
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1  Introduction 


Current  vibration-bcised  techniques  used  for  gearbox  fault  detection  rely  on  analysis  of 
signals  captured  under  fixed  load  and  fixed  speed  conditions.  While  helicopter  gearboxes 
operate  nominally  at  fixed  speed,  they  do  not  operate  under  fioced  load.  If  such  fault 
detection  techniques  are  to  be  implemented  as  part  of  a  self-contained  on-board  system, 
then  adaption  to,  or  accounting  for,  load  variation  is  a  necessary  requirement. 

Initial  investigations  of  gearbox  vibration  signals  indicate  that  load  variation  acts  to 
change  the  signal  characteristics  in  a  non-linear  manner.  This  implies  that  spectral  and 
statistical  properties  of  the  vibration  signal  vary  non-linearly  with  load.  Consequently, 
adaption  to  load  variation  cannot  be  achieved  by  simply  normalising  the  signal  with  respect 
to  the  load,  or  some  signal  feature  which  is  proportional  to  load. 

This  report  details  signal  processing  methods  which  are  able  to  account  for  load  vari¬ 
ation. 

One  of  the  difficulties  in  attempting  to  account  for  load  variation  in  the  vibration 
signal  is  that  it  is  not  possible  to  use  a  strictly  adaptive  method.  Such  methods  usually 
derive  an  error  signal  from  the  actual  vibration  signal,  and  a  mathematical  model  used 
to  predict  the  vibration  signal.  This  error  is  then  used  to  adjust  the  model  parameters  in 
order  to  reduce  the  error  power.  A  problem  arises  when  a  fault  occurs  in  the  system,  in 
which  case,  the  model  will  be  adjusted  not  only  to  account  for  the  variation  in  load,  but 
also  to  any  change  in  system  dynamics  due  to  the  fault.  The  aim,  then,  is  to  account  for 
load  variations  while  maintaining  the  ability  to  detect  faults. 


2  Signal  Model 

The  signal  processing  methods  used  herein  essentially  involve  modelling  the  signal  of 
interest,  and  adjusting  the  model  parameters  appropriately,  to  account  for  varying  load 
conditions.  The  underlying  model  may  be  either  (1)  linear  or  (2)  non-linear.  While  the 
theory  of  linear  models  is  well  established,  non-linear  models  have  found  prominence  in 
recent  years  [7,  9,  10,  11,  15,  24,  25]. 

There  are  many  types  of  non-linear  models  including  non-linear  networks,  bilinear 
models  and  Volterra  series.  Due  to  their  complexity,  and  the  computational  burden,  non¬ 
linear  models  should  only  be  used  if  absolutely  necessary.  Conversely,  linear  models  are 
much  easier  to  implement,  and  require  less  computational  effort.  In  this  report  only  linear 
models  will  be  considered,  however,  the  techniques  are  readily  extended  to  non-linear 
models. 

One  of  the  most  general  linear  models  is  the  autoregressive  moving  average  (ARMA) 
model.  Such  models  have  been  used  in  fields  as  diverse  as  economic  forecasting,  seismology, 
speech  processing,  array  processing  and  adaptive  filtering  [13,  14,  16,  19,  20,  21,  22]. 
ARMA  models  are  also  widely  used  in  control  theory. 
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2.1  General  Linear  Signal  Model 


The  signal  under  investigation  is  recorded  from  an  accelerometer  mounted  externally 
on  the  gearbox.  Since  the  measured  vibration  signal  is  stochastic  (as  opposed  to  determin¬ 
istic),  it  is  convenient  to  model  it  as  a  white  noise  process  passed  through  a  linear  filter 
with  a  rational  spectral  density  function  as  shown  in  figure  1.  Here  {x[n]  :  n  =  1,2,. . . .  N} 


White  noise 
input 


v[n] 


Vibration 

signal 


2:[n] 


Figure  1:  Vibration  signal  model. 

is  the  recorded  vibration  signal  sequence  and  {w[n]  :  n  =  1, 2, ....  iV}  is  a  Gaussian  dis¬ 
tributed  random  sequence  with  zero  mean  and  variance  (often  it  is  assumed  <7^  =  1  for 
convenience).  The  polynomials  A{z)  and  B{z)  represent  the  poles  and  zeros,  respectively, 
of  the  linear  system.  If 

A{z)  =  1  (1) 

the  model  consists  entirely  of  zeros  and  is  known  as  a  moving  average  (MA)  process. 
Alternatively,  if 

B{z)  =  1  (2) 

the  model  consists  entirely  of  poles,  and  is  known  as  an  autoregressive  (AR)  process.  Both 
A{z)  and  B{z)  are  assumed  to  be  polynomials  in  z  (the  forward  shift  operator)  of  finite 
order,  and  share  no  common  factors. 

It  is  generally  accepted  [4,  13,  14,  16]  that  AR  models  are  suited  to  processes  which 
exhibit  peaks  in  the  power  spectral  density  function,  and  MA  models  are  suited  to  processes 
which  exhibit  notches  in  the  power  spectral  density  function.  Although  the  AR  and  MA 
models  (individually)  are  able  to  represent  a  broad  range  of  signals  over  an  arbitrary 
large  frequency  range,  the  ARMA  model  is  able  to  achieve  this  with  a  smaller  number  of 
parameters. 

It  is  important  to  note  a  fundamental  difference  between  signal  processing  and  control 
applications;  viz.,  the  properties  of  the  signal  v[n\.  In  general  signal  processing  applica¬ 
tions,  access  to  v[n]  is  unavailable,  and  it  is  assumed  to  consist  of  a  Gaussian  distributed 
process  as  described  above.  In  control  applications,  however,  u[n]  is  often  a  deterministic 
function,  of  which  access  is  generally  available.  Access  to  v[n]  makes  the  identification  of 
the  model  B{z)/A{z)  simpler  than  if  it  is  not  available.  Additionally,  for  Gaussian  n[n], 
B{z)/A{z)  is  assumed  to  be  a  stable  minimum-phase  system  (poles  and  zeros  are  inside 
the  unit  circle). 


2.2  Fault  Signal  Model 

The  linear  signal  model  above  describes  the  process  under  normal  operation  (free  from 
any  fault).  The  fault  may  be  incorporated  by  adding  a  signal  which  describes  the  action 
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of  the  fault, 


x\n]  =  j/[n]  +  a  z[n] 


(3) 


where  j/[n]  is  the  signal  recorded  during  normal  conditions,  z[n]  is  the  signal  due  to  a  fault, 
and  a  =  0  if  normal  operation  or  a  =  1  if  operation  under  fault  conditions.  Figure  2  shows 
this  graphically.  The  sequence  {w[n]  :  n  =  1, 2, . . . ,  N}  is  Gaussian  distributed  with  zero 


Figure  2:  Fault  signal  model. 

mean  and  variance  Additionally,  {u[n]}  and  {w[n]}  are  assumed  to  be  independent. 
The  polynomials  C{z)  and  D{z)  describe  the  spectral  characteristics  of  the  fault  signal 
sequence. 


3  Accounting  for  Load  Variation 

The  simplest  method  of  accounting  for  load  variation  is  to  adjust  the  model  parameters 
with  respect  to  the  load.  This  is  illustrated  in  figure  3. 

Load 


t;[n] 


Figure  3:  Adaptive  parameter  scheme. 


This  requires  the  model  parameters  to  be  identified  as  a  function  of  the  load, 
model  may  then  be  written  as 


H{<:,z)  = 


MC,z) 


where  the  variable  C  is  introduced  to  represent  the  dependence  on  load. 


The 

(4) 


This  scheme  requires  either  the  parameter  variation  to  be  determined  as  a  function 
of  load,  or  more  simply,  parameters  to  be  identified  for  a  series  of  different  loads  C[^]  = 
Bmini  •  ■  •  )  Lmax  where  Lmin  and  Lmax  are  the  minimum  and  maximum  loads  of  interest. 
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It  will  be  necessary  to  determine  these  parameters  from  vibration  signals  which  are  known 
to  be  fault  free. 


4  Parameter  Estimation 


One  of  the  important  requirements  in  using  the  modelling  techniques  described  above 
is  estimating  the  model  parameters.  These  model  parameters  must  be  estimated  from  the 
gearbox  vibration  signal.  There  are  several  techniques  for  ARMA  model  estimation,  and 
in  the  following  four  such  methods  are  described. 

It  is  well  established  [4,  13,  14,  16]  that  the  estimation  of  the  AR  parameters  involves 
a  set  of  linear  equations,  whereas  the  estimation  of  the  MA  parameters  is  a  non-linear 
problem.  The  latter  of  these  is  the  more  difficult  to  solve. 

The  first  three  methods  presented  below  differ  in  the  way  they  estimate  the  MA  pa¬ 
rameters  of  the  model.  They  are  also  static  estimators  as  they  operate  on  a  single  data 
block.  The  fourth  technique  uses  a  recursive  least  squares  algorithm  to  estimate  the  AR 
and  MA  parameters  concurrently.  It  is  also  an  adaptive  algorithm  since  it  updates  the 
parameters  continuously. 


4.1  Kay  and  Marple  Methods 

The  books  by  Kay  and  Marple  [13,  16]  present  several  methods  of  ARMA  parameter 
estimation.  Although  the  general  motivation  is  spectral  estimation,  these  methods  are 
also  applicable  to  system  modelling  and  linear  prediction.  Two  methods  are  presented  in 
which  the  AR  and  MA  parameters  of  the  ARMA  model  are  estimated  separately.  These 
are  referred  to  as  the  modified  Yule- Walker  equations  and  the  least  squares  modified  Yule- 
Walker  equations. 

The  general  form  of  a  rational  polynomial  ARMA  model  is 

x[n]  =  —  ^  a[k]x[n  —  k]  +  '^  b[k]v[n  —  k],  (5) 

fc=i  fc=l 


where  v[n]  is  a  Gaussian  distributed  white  noise  process,  and  a[k]  :  A:  =  1, 2, . . .  ,p  are  the 
autoregressive  parameters  and  b[k]  :  k  =  1,2, . . .  ,q  are  the  moving  average  parameters  for 
and  ARMA(p,  q)  process.  Multiplying  equation  (5)  by  x[n  —  k]  and  taking  the  expected 
value  gives 


Rxx[k]  =  -  E  - 1]  +  E  bm^vik  - 1], 


(6) 


1=1 


1=1 


where  Rxx[k]  is  the  autocorrelation  of  the  output  sequence  x[n]  and 


RxuM  =  E  {u[n]a;[n  -  fc]}  =  E  {x[n]v[n  -f  A:]}  . 


(7) 
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Now  Rxv[k]  =  0  for  fc  >  0,  since  x[n]  =  Ya=-oo  output  of  a  causal  filter 

and  clearly  depends  on  {u[n],  v[n  —  1], . . .}  which  is  uncorrelated  with  v[n  +  k]  for  A:  >  0. 

p  q 

-  a[l]Rxx[k  -  ^]  +  X!  -1]  for  A;  =  0, 1, . . . ,  g 

i2xxW  =  |  (8) 

—  Y,  <^mRxx[k  —  /]  for  A;  >  g  +  1 

.  1=1 

The  above  relationship  shows  there  is  a  non-linear  relationship  between  the  autocorrelation 
function  and  the  ARMA  parameters.  This  is  due  to  the  b[l]Rxv[k  —  i]  term.  If  only 
the  equations  for  A:  >  5  are  used,  then  this  term  may  be  omitted.  This  may  then  be 
written  £is  the  set  of  linear  equations 

All  [5]  Axx[9~1]  •••  Ais  [9  —  P  f]  ^[f]  Aii[9  + 1] 

Axi[9  + 1]  Rxx[q]  •••  Aii[9”P  +  2]  0(2]  Rxx\ij  +  ^ 

Rxx[q+p-^  Rxx[q  +  p-2]  AxxM  J  [  “W  J  lRxx[q  +  p] 

These  are  known  as  the  modified  Yule- Walker  equations.  These  equations  involve  the 
actual  autocorrelation  function  which  is  not  generally  known.  A  reasonable  approach  is 
to  replace  the  theoretical  autocorrelation  function  by  an  either  the  biased  estimate 

1 

Rxx[k]  =  ji^Y  -  *],  (10) 

i—1 

or  the  unbiased  estimate 

N-k 

^xx[k]  =  Y  -  ?].  (11) 


4.1.1  Modified  Yule- Walker  Equations 

The  modified  Yule- Walker  equations  may  be  solved  by  using  an  extension  of  the  Levin¬ 
son  [13,  16]  algorithm  to  give  estimates  of  the  AR  parameters  {a[l],a[2], . . .  ,a[p]}.  The 
recursive  algorithm  is  initialised  by  setting 

fii  ~Axi[9  +  1] 


i,,[i]  = 


~Rxx\q  1] 

Rxx[q] 


Pi  =  (1  -  ai[l]6i[l])  Ra;2;[g], 

Then  recursion  for  A:  =  2, 3, ...  ,p, 

r,  1  _  ~Rxx[q  A:]  +  Yli=i  0‘k—ib]Rxx[q  -l-  A:  —  ^] 


ajtW  =  ak-i[i]  +  ak[k]bk-i[k  -  i]  i  =  1, 2, . . . ,  A:  -  1. 
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If  fc  =  p,  the  recursion  is  stopped,  otherwise 

~Rxx\(l  ~  ^k  —  \\J'\Rxx\(l  —  k  —  l\ 


bk[k]  = 


Pk-1 


OjfcW  =  -*]  i  =  1,2, . . . ,  k  -  1, 

Pk  =  (1  -  aA:W6fc[*])Pfc-l- 


(14) 


The  solution  is  then  a[k]  =  ap[k]  :  k  =  1,2,. . .  ,p,  and  the  error  variance  =  Pp. 

Following  the  estimation  of  the  AR  parameters,  the  data  are  filtered  with  the  estimate 
of  A(z), 

A{z)  =  a[l]z~^  +  a[2]z”^  +  . . .  +  a[p]z~'^,  (15) 


giving, 

y{^]  =  “  Z]  ~  (16) 

fc=i 

Durbin’s  method  [13,  16]  is  applied  to  this  filtered  data  to  estimate  the  MA  parameters 
{5[lj,  5[2], . . . ,  6[g]}  and  the  white  noise  variance  a^. 

Durbin’s  method  consists  of  the  following  two  steps, 


1.  Using  the  data,  {p[l],  y[2], . . . ,  yfAT]},  fit  a  large  order  AR  model  using  the  autocor¬ 
relation  method  (based  on  the  Levinson  algorithm).  For  an  AR  model  of  order  L, 
where  q  L  N,  the  white  noise  variance  estimator  is  given  by 

L 

=  JSx2;[0]  +  53  a.[k]Rxx[k].  (17) 

*=1 


2.  Using  the  AR  parameter  estimates  obtained  in  step  1  as  the  data  {a[l],  a[2], . . . ,  a[L]}, 
use  the  autocorrelation  method  with  an  order  q  to  find  {5[1],  5[2], . . . ,  b[q]}. 


4.1.2  Least  Squares  Modified  Yule- Walker  Equations 

To  reduce  the  variance  associated  with  the  modified  Yule- Walker  equation  estimator 
described  above,  it  is  possible  to  utilise  more  of  the  available  equations.  Since  for  an 
ARMA(p,  g)  process 


=  A:  >  g  + 1  (18) 

Z=1 

the  choice  of  only  using  p  equations  in  (9)  is  arbitrary.  The  motivation  for  extending 
this  number  of  equations  is  that  the  autocorrelation  function  contains  information  beyond 
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the  p  lags  used  in  equation  (9)  [4] .  Assuming  that  the  highest  lag  of  the  autocorrelation 
function  which  may  be  estimated  is  M,  consider  the  following  set  of  equations, 


Rxx[q  + 1] 
■Rxxi?  +  2] 

Rxx\q\  Rxx\.q  1] 

Rxx[q  T 1]  -RxxM 

■  Rxx[q-p  +  l] 

•  Rxx[q-p  +  2] 

r  «[i]  1 

n[2] 

Rxx[M] 

•  Rxx  [At  p] 

.  "W  . 

(19) 


or  in  matrix  form 


r  =  — Ra. 


(20) 


With  R  of  dimension  {M  —  q)  x  p,  and  assuming  M  —  q  >  p,  there  will  be  more  equa¬ 
tions  than  unknowns.  Also,  since  the  actual  autocorrelation  is  replaced  by  an  estimate, 
allowance  must  be  made  for  estimation  error.  If  equation  (20)  is  re-written  as 


f  =  — Ra  -t-  e. 


(21) 


then  a  least-squares  solution  for  the  AR  parameters  [13,  18]  is  given  by, 

a  =  -  (R^R)“^R^f. 


(22) 


The  estimator  is  asymptotically  unbiased  and  has  variance  which  usually  decreases  as 
M  —  q  (the  number  of  equations)  increases,  provided  N  ^  M  —  q. 

To  determine  the  MA  parameters  and  white  noise  variance  estimate,  Durbin’s  method 
may  be  used  (as  described  above),  again  after  filtering  with  the  estimate  of  A(z). 


4.2  Mayne  Firoozan  Method 

The  Mayne-Firoozan  method  of  ARMA  parameter  estimation  [17]  fits  a  high  order 
AR  model  to  the  data,  then  filters  the  data  with  the  inverse  of  this  model  as  a  means 
of  estimating  the  white  noise  process.  This  filtered  data  are  then  used  in  conjunction 
with  the  original  data  to  determine  an  ARMA  input/output  model  via  a  least  squares 
estimation. 

The  process  consists  of  5  steps, 

1.  Determine  A{z),  the  value  of  A{z)  which  minimises 

N  N 

X)  [^iz)y[k]f  =  [y[k]  +  aiy[k  -  1]  +  ...  -1-  aMy[k  -  M]f  ,  (23) 

fc=l  A:=l 

where  M  ^  p. 

2.  Determine  the  residual  sequence  e[k]  :  A:  =  1, 2, . . . ,  W  using 


e[k]  =  A{z)y[k]. 


(24) 
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3.  Determine  A^{z)  and  B^{z),  the  values  of  A{z)  and  B{z)  which  minimise 

N  N 

Y^[A{z)y[k]-B{z)e[k]]  =  ^  [(y[A:]  -  e[A:])  +  (ci, . . .  ,ap, 6i, . . . 6p) 

x(?/[A:-  l],...,y[A;-p]e[A:-  1], . . .  ,e[A;  - 

4.  Determine  the  filtered  sequences  {y}  and  {e}  using 

B{z)y[k\  =  y[k], 

B{z)e[k]  =  e[k]. 

5.  Determine  A‘^{z)  and  B^{z),  the  values  of  A{z)  and  B(z)  which  minimise 


N 


N 


'Y^[A{z)y[k]  -  B{z)e[k]]  =  ^  [(y[A;]  -  e[fc])  +  (oi, . . .  ,ap,  6i, . . .  6p) 


k=l 


k=l 


X -  1], . . . , y[A:  -  p]e[k  -  1], . . . ,  e[^-  -  p])] 


(25) 


(26) 


(27) 


Steps  1,  3  and  5  are  accomplished  using  a  linear  least  squares  estimate.  Steps  4  and  5 
are  added  to  improve  the  consistency  of  the  parameter  estimates. 


4.3  Stochastic  Recursive  Least  Squares 


The  recursive  least  squares  algorithm  is  used  widely  in  adaptive  control  [1,  2,  12]  and 
adaptive  filtering  [8],  as  a  method  of  continuously  updating  model  parameters  as  new 
data  are  received.  In  both  these  fields,  however,  access  is  available  to  the  system  (or 
filter)  inputs  and  outputs.  As  noted  previously,  in  vibration  analysis  there  is  generally  no 
input  signal  available,  thus  the  least  squares  algorithm  requires  a  slight  modification  [12] 
to  produce  a  stochastic  recursive  least  squares  algorithm.  The  algorithm  is  summarised 
below,  and  a  detailed  derivation  is  presented  in  Appendix  A. 

The  parameterised  model  is  written  in  the  form 

y[n]  =  ipi [n]9i  +  ip2[n]62  +  . . .  +  ipM[n]9M,  (28) 


where  the  model  parameters  are, 

9  =  [aia2...apbib2...hp],  (29) 

ie.,  M  =  2p.  The  system  outputs,  y[n],  and  estimated  inputs,  i;[n]  are  augmented  to  form 
the  array, 

(p^[n]  =  [y[n  —  l]y[n  —  2] . . .  y[n  —  p]v[n  —  l]u[n  -  2] . . .  v[n  —  p]] .  (30) 

The  array  of  outputs  used  in  the  above  model  up  to  time  n  is  then 


y[n]  = 


y[i] 

y[2] 

yW 


(31) 
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and  the  model  error  is  defined  as, 


e[n]  =  y[n]  -  y[n\ 

=  y[n]  -  (P^[n]d, 


(32) 


where  y[n]  us  the  estimate  of  y[n]  based  on  the  estimated  parameters  6.  The  array  of 
errors  may  be  defined  as 


E[n]  = 


e[l] 

e[2] 

efn] 


=  Y[n]  -  (f}[n]§, 


(33) 


where, 


(j)[n]  = 

1 

1 _ 

with  ip[n\  = 

ipi[n] 

'P2[n] 

L  J 

.  <^p[n]  _ 

(34) 


Or  in  terms  of  the  outputs  and  estimated  inputs, 

4>[n]  = 


-3/[0] 

-2/[l] 


-y[o] 


-y[l-p]  ^;[0]  v[-l] 

-y[2-p]  ?;[1]  t;[0] 


?;[!  -p] 

v[2  -p] 
v[n  —  p] 


-y[n  —  1]  —y[n  —  2]  ...  —y[n  —  p\  v\n  —  1]  v{n  —  2] 

For  simplicity  in  the  algorithm,  define  the  matrix 

P[n]  =  4>[n^  . 

The  least  squares  recursive  algorithm  for  estimating  the  parameter  vector  9  is 
1.  Determine  K[n\  (refer  to  Appendix  A  for  details), 

K[n\  =  P[n  —  l]^[n]  +  ip^[n]P[n  -  l](^[n]j  , 


(35) 


(36) 


(37) 


where  I  is  the  identity  matrix.  For  a  single  input  single  output  system,  this  reduces 
to 

P{n  —  l]¥?[n] 


1  +  (p'^\n\P[n  —  l](/3[n] ' 

The  matrix  K\p\  is  often  referred  to  as  the  Kalman  gain  matrix  [8]. 
2.  Update  the  estimated  parameter  vector  6 

0[n]  =  6[n  —  1]  +  K[n]  |^y[n]  —  (p^[n\9[n  —  l]j  . 


(38) 


(39) 
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3.  Update  the  matrix  P[n] 

P[n]  =  [I  -  K[n]ip[n]]  P[n  -  1]. 
The  procedure  is  initialised  by  setting 


P[0]  = 


a  0  ...  0 

0  a  ...  0 

:  :  0 
0  0  ...  a 


=  al, 


(40) 


(41) 


for  some  large  a. 

The  above  algorithm  assumes  knowledge  of  the  input  sequence,  which  is  unavailable  for 
a  stochastic  model.  However,  since  the  prediction  error  for  a  stochastic  model  is  analogous 
to  the  Gaussian  distributed  input  sequence,  it  may  be  used  as  an  estimate  of  the  input 
sequence, 

v[n]  =  e[n]  =  y[n]  -  (p^[n]9[n  -  1].  (42) 

To  initialise,  the  error  sequence  is  assumed  zero  for  n  <  0  ie.,  v[n]  =  0.  Additionally,  the 
initial  parameter  vector  is  also  set  to  zero,  0[O]  =  0. 


5  Fault  Detection  Condition  Indices 

In  order  to  determine  the  presence  or  absence  of  a  fault,  some  signal  feature  or  condition 
index  must  be  determined.  These  are  generally  a  numerical  value  derived  from,  say,  the 
statistics  of  the  vibration  signal.  Empirical  analysis  and  experience  allow  ranges  of  these 
values  for  which  acceptable  and  unacceptable  operation  may  be  defined. 


5.1  Current  Methods 

At  AMRL,  the  use  of  synchronous  signal  averaging  has  been  fundamental  to  the  ability 
to  detect  faults  on  individual  gears  within  a  gearbox.  The  approach  used  is  to  average 
segments  of  the  vibration  signal  which  correspond  to  one  shaft  rotation.  Such  ensemble 
averaging  allows  signal  frequencies  which  are  synchronous  with  the  shaft  period  for  that 
gear  to  be  isolated  from  all  other  signals. 

The  signal  average  contains  a  large  amount  of  information  which  is  redundant  to  the 
aim  of  fault  detection.  This  is  essentially  contained  in  the  gear  meshing  frequency  and  its 
harmonics,  as  well  as  a  smaller  amount  at  the  fundamental  shaft  rotation  frequency,  and 
its  modulation  with  the  meshing  harmonics.  The  latter  may  be  due  to  imbalances  in  shaft 
rotation  (perhaps  caused  by  shaft  misalignment).  However,  this  once-per-rev  component 
may  be  a  vital  factor  once  localised  faults  propagate  sufficiently  as  shown  by  Forrester  [6] . 
For  example,  a  defect  on  a  single  gear  tooth  will  have  a  prominent  once-per-rev  component 
each  time  it  comes  into  mesh,  provided  the  effect  of  this  fault  exceeds  that  due  to  shaft 
imbalances.  Thus  there  is  a  trade-off  in  removing  the  once-per-rev  component. 
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The  information  is  redundant  in  terms  of  fault  detection  since  it  is  relevant  to  the  gear 
meshing  only,  and  contains  no  information  relevant  to  the  fault.  Hence  it  may  be  removed 
as  a  means  of  enhancing  sensitivity  to  faults. 

Removal  of  the  tooth  meshing  harmonics,  as  well  as  the  once-per-rev  and  possibly 
twice-per-rev  modulation  of  these,  produces  a  residual  signal  average.  Clearly  this  residual 
signal  average  would  be  more  useful  for  fault  detection  analysis,  and  some  fault  indices 
used  at  AMRL  are  based  on  this  premise. 

Current  condition  indices  used  at  AMRL  [3,  5]  have  been  found  by  experience  to  be 
useful  indicators  of  various  types  of  gearbox  faults.  The  indices  are  of  two  basic  types:  (1) 
general  and  (2)  local.  General  indices  are  sensitive  to  general  faults,  such  as  wear  on  all 
teeth;  whereas  local  indices  are  sensitive  to  localised  perturbations,  such  as  changes  in  the 
stiffness  of  a  single  tooth  on  a  gear.  The  indices  used  at  AMRL  are  summarised  below. 

B1  Amplitude  of  the  shaft  frequency. 

A1  Amplitude  of  2  x  shaft  frequency. 

G1  RMS  of  the  signal  average. 

G2  RMS  of  the  residual  signal  average  divided  by  the  RMS  of  the  signal  average. 

G3  The  peak  level  of  the  envelope  of  the  narrow  band  around  a  mesh  harmonic  divided 
by  the  RMS  of  the  mesh  harmonic. 

LI  Kurtosis  of  the  residual  signal. 

L2  Kurtosis  of  the  envelope  of  the  narrow  band  around  a  mesh  harmonic. 

L3  Kurtosis  of  the  instantaneous  frequency  of  the  modulated  signal. 

5.2  Methods  for  Diagnosis 

As  described  above  the  use  of  synchronous  signal  averaging  and  removal  of  gear  mesh 
harmonics  allows  the  removal  of  a  large  amount  of  redundant  information  (at  least  in 
terms  of  fault  detection).  The  signal  averaging  process  enables  frequencies  relevant  to 
the  particular  gear  of  interest  to  be  extracted  from  those  due  to  other  gears  (provided,  of 
course,  other  gears  do  not  have  rotation  frequencies  and  harmonics  common  to  that  gear) . 
Further  to  this,  the  residual  signal  average  removes  the  gear  meshing  harmonics  which  are 
of  little  interest  to  fault  detection. 

With  this  in  mind,  and  the  use  of  ARMA  modelling  as  a  means  of  accounting  for 
load  variation,  the  residual  signal  obtained  from  the  model  appears  intuitively  to  offer  the 
most  promise  as  a  means  of  fault  detection.  The  residual  signal,  in  this  case,  is  just  the 
difference  between  the  measured  vibration  signal,  and  the  one-step-ahead  prediction  based 
on  the  estimated  ARMA  model, 

e[n]  =  a:[n]  —  :c[n] 

=  a:[n]  -1-  ^  a[A:]a:[n  —  fc]  —  S[fc]u[n  —  k] 

k=l  k=l 
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It  should  be  stressed  that  the  residual  signal  here  is  not  the  same  as  the  residual  signal 
average^ . 

Since  the  ARMA  model  attempts  to  model  the  dynamics  of  the  process,  the  residual 
signal  is  ideally  a  Gaussian  distributed  random  process  of  variance  ie.,  it  is  equivalent 
to  the  input  of  the  ARMA  model.  To  clarify  this,  the  ARMA  model  is  defined  as  a 
Gaussian  process  passed  through  a  linear  filter  of  the  form  B{z)IA{z).  Using  this  model 
to  predict  the  output,  the  error  must  be  equivalent  to  the  noise  input  since  this  is  the  part 
of  the  signal  which  cannot  be  predicted. 

Once  a  fault  is  introduced,  the  structure  of  the  process  model  changes.  Thus  using 
a  predictor  based  on  a  fault-free  process  will  change  the  residual  signal  from  Gaussian 
noise  to  coloured  (filtered)  noise  since  there  will  be  un-modelled  dynamics  due  to  the 
introduction  of  the  fault. 

From  a  dififerent  perspective,  but  conceptually  the  same,  the  residual  signal  is  akin  to 
the  process  of  synchronous  signal  averaging  and  removal  of  the  gear  mesh  harmonics  as 
described  above.  The  most  notable  difference  is  that  the  latter  uses  ensemble  averaging 
whereas  the  former  does  not.  The  residual  signal  consists  of  the  original  signal  minus 
the  identified  process  dynamics.  Hence  removal  of  the  redundant  information  is  similarly 
achieved. 

The  residual  signal,  then,  provides  the  most  obvious  means  of  detecting  the  presence 
of  a  fault.  Statistical  properties  and  the  power  spectral  density  offer  the  best  means  of 
extracting  information,  from  this  signal,  which  may  be  used  to  detect  faults. 

Coherent  signal  averaging  offers  another  advantage  hitherto  unmentioned.  That  is  it 
provides  a  convenient  description  of  the  signal  from  a  particular  gear.  In  a  crude  sense, 
ensemble  averaging  compresses  the  information  recorded  over  many  minutes  into  a  short 
data  record,  which  describes  the  rotation  of  a  single  gear,  and  this  is  primarily  what  we 
are  interested  in.  Application  of  synchronous  signal  averaging  to  the  residual  signal  from 
the  ARMA  model  will  then  enable  extraction  of  information  pertinent  to  the  rotation  of  a 
single  gear,  particularly  if  localised  faults  are  present,  since  (ideally)  all  other  information 
has  been  removed.  In  this  case,  current  fault  indices  based  on  the  residual  signal  average 
may  be  use  on  the  signal  average  derived  from  the  ARMA  model  residual  signal. 

The  use  of  ARMA  modelling,  then,  enables  load  variation  to  be  accounted  for,  and  the 
residual  of  the  model  provides  the  means  of  detecting  faults. 


6  Simulations 

In  this  section  a  simple  ARMA  time-series  is  generated  to  simulate  a  linear  process. 
To  simulate  varying  load  conditions,  the  poles  and  zeros  of  the  model  will  be  altered  in 
steps.  To  simulate  the  action  of  a  fault,  a  second  ARMA  time-series  is  generated,  and 

^The  residual  signal  average  is  the  ensemble  average  of  the  vibration  signal  with  the  meshing  harmonics 
removed.  The  residual  signal  in  this  case  is  the  output  of  a  prediction  error  filter,  and  no  ensemble  averaging 
is  applied. 
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added  to  the  first  for  a  short  period  of  time.  Since  the  dynamics  of  the  fault  are  also  likely 
to  change  with  varying  load,  the  poles  and  zeros  of  the  fault  signal  are  also  varied  with 

load. 

The  following  simulations  are  intended  only  as  an  illustration  of  how  load  variation  may 
be  accounted  for,  and  how  fault  detection  is  possible.  Further  work  using  data  recorded 
from  the  AMRL  gear  test  rig  will  appear  in  a  subsequent  report. 


6.1  General  Time  Series  Simulation 


To  generate  ARMA  time-series  examples  with  a  known  transfer  function,  it  is  necessary 
to  separate  the  AR  and  MA  sections  of  the  model.  Consider  the  ARMA  time-series  model 
of  equation  (5)  ^ 

x[n\  —  —  ^  a[k]x[n  —  /s]  + 

it=l  k=l 

where  :  A:  =  1,2, . . .  ,p}  is  a  Gaussian  distributed  white  noise  process.  Since  the 

AR  section  requires  knowledge  of  x[n]  for  n  <  0  (and  this  is  clearly  not  available)  it  is 
necessary  to  determine  appropriate  initial  conditions.  This  is  achieved  by  setting  [13] 


y[n]  = 


^^^ov[0] 

k 

-Y^an[l]y[n-l]  + ^/^v{n] 
-  1=1 


for  rz  —  0 

for  n  =  1, 2, . . .  ,p  —  1 


(45) 


where  Oi[i]  is  the  zth  prediction  coefficient  for  the  jth  order  model,  and  pj  is  the  prediction 
error  power  for  the  jth  order  model.  This  requires  an  additional  p  1  sets  of  prediction 
coefficients  (for  the  predictors  of  order  1, 2, ...  ,p  -  1)  which  may  be  derived  from  the  pth 
order  AR  model  by  [13] 


1  [*] 

Pk-l 


ak[i\  -  ak[k]ak[k  -  i] 


1  -  |afc[A:]p 


Pk 


\ak[kW 


for  k 


forz  =  l,2,...,A:  — 1  A:  =  p  — l,p  — 2,...,1 

p  -  l,p  -  2, . . . ,  1. 


(46) 


With  the  above  initial  conditions  set,  the  remaining  data  for  the  AR  model  may  be 
generated  using 


p 

y[n]  = -'^a[n]y[n-k]  + y/^v[n]  n  =  p,p -j- 1, . . . ,  iV -f  g,  (47) 

fc=l 

where  N  is  the  total  number  of  time-series  samples  required. 

Filtering  the  AR  time-series  {y[n]  :  n  =  1,2,  with  a  finite  impulse  response 

(FIR)  filter  consisting  of  coefficients  equal  to  the  MA  coefficients  will  then  produce  the 
required  ARMA  time-series, 


k=l 


n  =  q,q  +  l,-..,N  +  q. 


(48) 
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Note  that  a  delay  of  q  samples  is  required  to  ensure  the  filter  is  initialised. 

This  ARM  A  time-series  is  completely  described  by  its  ARM  A  parameters  {a[fc]  :  k  = 
1, 2, . . .  ,p}  and  {6[A:]  :k  =  1, 2, . . . ,  q}  as  well  as  the  prediction  error  power  pp. 


6.2  System  Model  and  Fault  Simulation 


Since  the  fundamental  frequency  components  of  the  vibration  signal  due  to  gear  mesh¬ 
ing  are  directly  related  to  the  shaft  rotation  frequency,  then  load  variation  will  not  change 
these.  Hence,  the  argument  (angle)  of  the  model  poles  and  zeros  in  the  .2-domain  will 
remain  the  same.  It  is,  however,  possible  that  the  damping  of  the  poles  and  zeros  (distance 
from  the  origin)  may  change  with  load. 

The  simple  system  model,  for  illustrating  how  load  variation  may  be  accounted  for, 
consists  of  two  pairs  of  complex  conjugate  poles,  and  a  double  real  zero,  as  shown  in 
figures  4  and  5.  Note  that  the  two  figures  show  the  initial  and  final  load  positions.  Also, 
the  double  zero  appears  only  once  on  the  real  axis.  To  simulate  load  variation  the  poles 
and  zeros  are  moved  from  the  initial  position  (figure  4)  to  the  final  position  (figure  5)  in 
five  steps. 


Figure  4:  Initial  system  poles  and  zeros. 


Figure  5:  Final  system  poles  and  zeros. 


Previous  analysis  of  vibration  signal  from  helicopter  transmission  systems  with  dam¬ 
aged  teeth  [5]  has  shown  that  the  fault  appears,  predominantly,  at  the  second  and  higher 
order  harmonics  of  the  tooth-meshing  frequency.  To  exaggerate  this,  for  illustration  pur¬ 
poses,  the  fault  signal  is  simulated  as  a  pair  of  complex  conjugate  poles,  and  a  single  zero 
in  the  left  hand  plane  of  the  .2-domain.  These  are  shown  in  figures  6  and  7.  Again  load 
variation  is  simulated  by  moving  the  poles  and  zero  away  from  the  origin  towards  the  unit 
circle  in  five  steps.  In  addition,  a  low  signal-to-noise  ratio  is  achieved  by  setting  the  white 
noise  variance  used  to  generate  the  fault  time-series  to  one  tenth  that  of  the  system  signal. 

The  spectra  of  the  system  and  fault  processes  are  shown  in  figures  8-11.  Note  that  the 
magnitude  of  the  system  spectra  are  significantly  greater  than  that  of  the  fault  spectra. 
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Figure  6:  Initial  fault  poles  and  zero.  Figure  7;  Final  fault  poles  and  zero. 


These  spectra  are  determined  from  the  actual  ARMA  parameters  using, 


\PARMA[f]\  - 


a^EUb[k]e-P-P 

1  +  ELi  a[k]e-P^f^  ■ 


(49) 


Figure  8:  Initial  system  ARMA  spectra.  Figure  9:  Final  system  ARMA  spectra. 


Figures  12  and  13  show  30,000  samples  of  the  system  ARMA  simulation  to  which  has 
been  added  10,000  samples  of  the  fault  signal  in  the  sample  range  10,000-20,000.  From 
these  figures  the  location  of  the  fault  is  not  apparent  in  either  of  the  simulations. 

Since  the  kurtosis  provides  a  measure  of  localised  faults,  a  series  of  impulses  with  a 
magnitude  of  10.0  is  added  to  the  fault  signal.  The  impulses  are  spaced  100  samples  apart. 
Time-series  simulations  for  the  initial  and  final  load  conditions  are  shown  in  figures  14  and 
15.  Again  the  location  of  the  fault  is  not  distinguishable,  despite  the  additional  impulses 
with  a  magnitude  of  10.0. 
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Figure  10:  Initial  fault  ARM  A  spectra. 
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Figure  11:  Final  fault  ARMA  spectra. 
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Figure  12:  Initial  system  and  fault  signal  Figure  13:  Final  system  and  fault  signal 
time-series.  time-series. 
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Figure  14:  Initial  system  and  fault  signal  Figure  15:  Final  system  and  fault  signal 
time-series  with  impulses.  time-series  with  impulses. 
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6.3  Fault  Detection 

To  test  both  general  and  local  fault  detection  indices,  the  variance  (equivalent  to  the 
RMS)  and  kurtosis  of  the  residual  signal  are  to  be  determined.  To  generate  the  residual 
signal,  a  prediction  error  filter  is  formed  using  the  ARMA  model.  A  diagram  is  shown  in 
figure  16. 


y[n] 


Figure  16:  Prediction  error  filter  block  diagram. 

For  the  data  generated  above  in  Section  6.2  an  ARMA(2,2)  model  is  estimated,  and 
a  prediction  error  filter,  as  shown  in  figure  16,  is  formed  using  the  parameter  estimates. 
Figures  17  and  18  show  the  residual  signals  determined  from  the  initial  and  final  fault 
signal  time-series  shown  in  figures  12  and  13.  In  figure  17  a  small  increase  in  signal 
magnitude  can  be  seen  from  10,000-20,000  samples  (the  location  of  the  fault),  however, 
the  increase  is  far  more  evident  in  figure  18.  This  suggests  that  some  measure  of  the  signal 
magnitude  will  he  suitable  for  detecting  a  fault.  The  reason  that  the  increase  is  greater 
for  the  final  fault  signal  is  that  the  poles  and  zeros  of  the  fault  signal  model  are  closer 
to  the  unit  circle.  This  has  the  effect  of  creating  a  sharper  notch  in  the  power  spectral 
density,  and  increasing  the  magnitude  of  the  fault  signal  (compare  figures  10  and  11,  and 
note  the  differences  in  the  magnitude  scales). 

Applying  the  same  prediction  error  filter  to  the  time-series  with  impulses  added  to  the 
fault  signals  give  the  results  shown  in  figures  19  and  20.  In  this  case  both  the  initial  and 
final  fault  signal  residual  signals  show  a  significant  increase  in  magnitude  at  the  location 
of  the  fault.  The  greater  discrimination  here  is  essentially  due  to  the  impulses,  which  axe 
not  accounted  for  by  the  ARMA  model  (since  they  are  part  of  the  fault  signal). 

Figures  21-25  show  the  variance  of  the  residual  signal  for  each  of  the  fault  signal  time- 
series  as  described  above.  The  variance  was  determined  by  averaging  over  1,000  samples 
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Figure  17:  Residual  signal  from  initial  fault  Figure  18:  Residual  signal  from  final  fault 
signal  time-series.  signal  time-series. 
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Figure  19:  Residual  signal  from  initial  fault  Figure  20:  Residual  signal  from  final  fault 
signal  time-series  with  impulses.  signal  time-series  with  impulses. 
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using  a  sliding  block  approach.  Here,  the  variance  is  determined  over  the  block  of  1,000 
samples.  Then  as  each  new  sample  is  added  to  the  start  of  the  block,  one  sample  is 
discarded  from  the  end,  and  the  variance  is  determined  for  the  new  data  block.  This  gives 
the  variance  over  the  averaging  time  (of  1,000  samples)  as  a  function  of  samples  (or  time). 

Note  that  despite  there  being  little  discrimination  between  normal  and  fault  conditions 
in  the  residual  signal  of  figure  17  (for  the  initial  fault  signal),  figure  21  shows  a  clear 
distinction  from  10,000-20,000  samples.  This  is  a  consequence  of  the  averaging  process 
used  in  determining  the  variance. 
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Figure  21:  Variance  of  residual  signal  from 
initial  fault  signal  time-series. 


Figure  22:  Variance  of  residual  signal  from 
second  fault  signal  time-series. 
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Figure  23:  Variance  of  residual  signal  from  Figure  2f:  Variance  of  residual  signal  from 
third  fault  signal  time-series.  fourth  fault  signal  time-series. 


Figures  26-30  show  the  variance  determined  for  the  residual  of  the  fault  signal  time- 
series  with  impulses  added.  It  is  evident  that  the  addition  of  impulses  allows  a  greater 
discrimination  between  normal  and  fault  conditions.  This  should  be  understood  by  com¬ 
paring  the  residual  signals  of  figures  17  and  19,  which  show  that  the  impulses  lead  to  a 
significant  increase  in  the  magnitude  of  the  residual  signal  when  the  fault  is  present. 

Figures  31-35  show  the  kurtosis  of  the  residual  signal  for  each  of  the  fault  signal  time- 
series  described  previously.  A  similar  sliding  block  approach,  as  described  previously,  is 
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Figure  25:  Variance  of  residual  signal  from  Figure  26:  Variance  of  residual  signal  from 
final  fault  signal  time-series.  initial  fault  signal  time-series  with  impulses. 
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Figure  27:  Variance  of  residual  signal  from  Figure  28:  Variance  of  residual  signal  from 
second  fault  signal  time-series  with  impulses,  third  fault  signal  tirne-series  with  impulses. 
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Figure  29:  Variance  of  residual  signal  from  Figure  30:  Variance  of  residual  signal  from 
fourth  fault  signal  time-series  with  impulses,  final  fault  signal  time-series  with  impulses. 
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also  used  with  averaging  over  1,000  samples.  In  each  of  these  figures  it  is  not  possible  to 
determine  the  location  of  the  fault,  since  the  kurtosis  appears  to  be  similar  throughout 
each  record.  There  is,  however,  a  distinct  increase  in  the  kurtosis  at  the  beginning  and 
end  of  the  fault  section  in  figure  35.  This  can  be  attributed  to  the  transition  from  normal 
to  fault  condition,  and  thence  from  fault  to  normal  condition.  This  feature  implies  that 
the  kurtosis  is  able  to  detect  transient  phenomena  in  the  residual  signal,  provided  it  is  of 
significant  magnitude.  The  transition  in  the  residual  signal  can  be  seen  clearly  in  figmre  18, 
where  the  fault  signal  firom  10,000-20,000  samples  has  a  much  greater  magnitude  than  the 
normal  signal  either  side  of  this. 


2.6-1  .  .  .  •  •  • 

0  5000  10000  15000  20000  25000  30000 

Samples 

Figure  31:  Kurtosis  of  residual  signal  from 
initial  fault  signal  time-series. 
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Figure  32:  Kurtosis  of  residual  signal  from 
second  fault  signal  time-series. 
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Figure  33:  Kurtosis  of  residual  signal  from  Figure  34:  Kurtosis  of  residual  signal  from 
third  fault  signal  time-series.  fourth  fault  signal  time-series. 


Figures  36-40  show  the  kurtosis  determined  for  the  residual  of  the  fault  signal  time- 
series  with  impulses  added.  In  contrast  to  figures  31—35  the  location  of  the  fault  can 
be  seen  clearly  in  the  sample  range  10,000-20,000.  Also  noticeable  are  sharp  peaks  at 
the  beginning  and  end  of  the  fault  section.  This  shows  that  the  kurtosis  is  suitable  for 
determining  faults  which  are  impulsive,  as  well  as  being  able  to  detect  the  transition  as 
the  fault  signal  is  added  to  the  normal  signal. 
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Figure  35:  Kurtosis  of  residual  signal  from  Figure  36:  Kurtosis  of  residual  signal  from 
final  fault  signal  time-series.  initial  fault  signal  time-series  with  impulses. 


Figure  37:  Kurtosis  of  residual  signal  from  Figure  38:  Kurtosis  of  residual  signal  from 
second  fault  signal  time-series  with  impulses,  third  fault  signal  time-series  with  impulses. 


Figure  39:  Kurtosis  of  residual  signal  from  Figure  40:  Kurtosis  of  residual  signal  from 
fourth  fault  signal  time-series  with  impulses,  final  fault  signal  time-series  with  impulses. 
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7  Further  Requirements 

This  report  is  intended  to  review  signal  processing  methods  which  may  be  suitable  for 
accounting  for  variable  load  conditions  when  analysing  gearbox  vibration  signals  for  fault 
detection.  The  application  of  these  methods  requires  much  further  study. 

Before  applying  these  methods,  however,  if  is  necessary  to  make  a  comprehensive  study 
of  vibration  signals  under  varying  load  conditions.  It  will  be  necessary  to  achieve  this  using 
signals  which  axe  known  to  be  fault-free.  It  is  important  to  gain  as  much  information  as 
possible  about  “normal”  condition  signals  particularly  if  little,  or  no,  information  is  to 
be  assumed  regarding  the  fault.  The  main  reason  for  this  is  that  it  is  not  possible  to 
characterise  all  types  of  faults,  since  there  are  perhaps  many  faults  for  which  we  have  no 
knowledge  of  [23]. 

The  methods  described  will  have  to  be  applied  to  actual  vibration  signals  to  determine 
the  suitability  of  the  types  of  models  presented.  Following  this,  it  will  be  necessary  to 
determine  how  effective  they  are  in  detecting  faults.  In  particular,  the  time  between 
detection  and  catastrophic  failure  will  be  a  most  important  indicator  here. 

Previously  it  was  suggested  that  model  parameters  will  need  to  be  determined  at 
several  different  load  settings.  Subsequently,  the  analysis  for  fault  detection  will  also  have 
to  be  achieved  at  different  load  settings.  However,  depending  on  the  variation  of  model 
parameters,  it  may  be  possible  to  interpolate  model  parameters  between  these  fixed  load 
points  to  achieve  a  more  practical  system,  able  to  be  used  at  any  load. 

An  obvious  extension  would  be  to  apply  synchronous  signal  averaging  to  the  ARMA 
residual  signal.  This  has  two  possible  benefits  namely:  (1)  it  will  be  possible  to  determine 
whether  the  residual  signal  contains  periodic  components  at  the  meshing  frequencies  - 
which  ideally  should  have  been  removed,  and  (2)  as  a  means  of  data  reduction.  This  latter 
point  will  be  of  benefit  if  on-board  hardware  is  to  be  used,  in  which  case  data  storage  will 
be  at  a  premium. 


8  Conclusions 

This  report  introduced  linear  modelling  techniques  as  a  means  of  accounting  for  vary¬ 
ing  load  conditions  experienced  in  helicopter  transmission  systems.  Several  methods  of 
ARMA  parameter  estimation  were  outlined  which  enable  modelling  of  the  vibration  sig¬ 
nals  recorded  from  transmission  systems.  With  the  ultimate  aim  being  the  detection  of 
faults  under  varying  load  conditions,  it  was  suggested  that  a  prediction  error  filter  be  used 
to  generate  a  residual  signal.  Subsequent  analysis  of  this  residual  signal  will  then  enable 
detection  of  faults. 

A  comparison  was  made  between  the  residual  signal  derived  from  the  prediction  error 
filter,  and  the  residual  signal  average  in  common  use  by  AMRL.  It  was  proposed  that  the 
former  is  able  to  remove  frequency  components  which  are  inherent  in  the  vibration  signal 
(including  those  synchronous  with  the  shaft  rotation  of  the  gear  of  concern,  and  other  gears 
in  the  gearbox),  whereas,  synchronous  signal  averaging  is  incapable  of  removing  any  signal 
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component  which  is  synchronous  with  the  shaft  rotation.  The  removal  of  the  meshing 
harmonics  and  once-per-rev  modulations  is  a  means  of  removing  information  redundant 
to  fault  detection,  from  the  signal  average.  However,  as  noted  previously,  this  may  also 
remove  information  specific  to  the  fault.  For  example,  the  once-per-rev  component  will  be 
evident  if  a  localised  gear  fault  (such  as  a  tooth  crack)  is  present,  in  addition  to  a  similar 
component  due  to  shaft  imbalances.  The  residual  signal  derived  from  the  prediction  error 
filter  does  not  suffer  this  disadvantage  since  the  model  is  estimated  from  a  vibration  signal 
which  is  known  to  be  fault  free.  Thus  the  residual  signal  should  only  contain  information 
which  is  relevant  to  the  fatdt. 

Simple  simulations  were  used  to  illustrate  the  ability  of  the  technique  to  detect  a  fault 
under  varying  load  conditions.  A  simulated  fault  signal  was  added  to  a  simulated  vibration 
signal  for  a  short  period.  It  was  found  that  the  residual  signal  derived  firom  the  prediction 
error  filter  had  an  increase  in  power  during  the  time  the  fault  signal  was  applied.  Two 
indices  were  determined  from  the  residual  signal;  the  (1)  variance  and  (2)  kurtosis.  It  was 
found  that  the  variance  was  sensitive  to  changes  in  signal  magnitude,  and  the  kurtosis  to 
transients  in  the  residual  signal. 

It  was  suggested  that  synchronous  signal  averaging  of  the  residual  signal  from  the 
prediction  error  filter  would  be  a  useful  means  of  data  reduction,  as  well  as  a  check  on  the 
effectiveness  of  the  methods  in  removing  information  redundant  to  fault  detection  (such 
as  the  gear  meshing  harmonics). 

It  was  noted  that  further  work  is  required  in  analysing  gear  rig  vibration  signals  to 
determine  the  applicability  of  the  modelling  techniques. 
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Appendix  A 

Stochastic  Recursive  Least  Squares  Derivation 


This  appendix  gives  a  derivation  of  the  recursive  least  squares  algorithm  adapted  for 
stochastic  signals.  It  is  drawn  from  several  sources  [1,  2,  12,  8]  in  an  attempt  to  provide 
a  more  concise  and  simplified  derivation. 

The  parameterised  model  is  written  in  the  form 

y[n]  =  ^p\[n]6i  +  +  •  •  •  +  (Al) 

where  the  model  parameters  are, 

6  =  [0102  . . .  Op6i62  ■  ■  •  ,  (-^2) 

ie.,  M  —  2p.  The  system  outputs,  y[n],  and  estimated  inputs,  v[n]  are  augmented  to  form 
the  array, 

<^[n]  =  [y[n  —  l]y[n  —  2] . . .  y[n  —  p]v[n  —  l]n[n  —  2] . . .  v[n  —  p]] .  (A3) 


The  array  of  outputs  used  in  the  above  model  up  to  time  n  is  then 

Y[n]  = 


y[i] 

y[2] 


and  the  model  error  is  defined  as. 


e[n]  =  y[n]  -  y[n] 

=  y[n]  -  p^[n]6, 


(A4) 


(AS) 


where  y[n]  us  the  estimate  of  y\n\  based  on  the  estimated  parameters  6.  The  array  of 
errors  may  be  defined  as 


where. 


E[n]  = 


e[l] 

e[2] 

e[nl 


=  y[n]  -  <i)[n]e, 


4{n]  = 

r  1 

<p^[2] 

with  (p[n]  = 

■  (pi[n]  ■ 
ip2[n] 

.  V’pW  . 

(A6) 


(A7) 
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Or  in  terms  of  the  outputs  and  estimated  inputs, 


0[n]  = 


-2/[0] 

-2/[l] 


-2/[0] 


-y[l-p]  u[0] 

-y[2-p]  v[l]  7;[0] 


z;[l  -  p] 
v  [2  -  p] 


(A8) 


-y[n  -  1]  -y[n  -  2] 


— y[n  —  p\  v[n  —  1]  v[n  —  2] 


v[n  —  p] 


The  least  squares  solution  requires  the  sum  of  the  errors  squared  to  be  minimised.  The 
sum  of  the  errors  squared  is  given  by  (note  that  the  time  index  is  dropped  for  convenience) 


V 


[yr  _  . 

yTy  _  yT^0 


[y-4 

—  (fP'Y  4-  d'^(f^4>6. 


To  minimise  V  with  respect  to  6,  differentiate  w.r.t.  6  and  set  to  zero, 

4.  [y^y  e'^cfY  + 

=  0-<f'Y-  fY  +  2f(t>9 
=  -2(jP'Y  +  2(jP^(j)9  =  0. 


(A9) 


(AlO) 


Now  the  second  derivative  (5/50^  d/d9)V  =  2(l>^(p  >  0,  thus  V  is  minimised.  The 
condition  for  a  least  squares  estimate  of  9  is  then 


2(t>^9  =  20^y 

9  =  (j>^Y. 


(All) 


A  recursive  algorithm  for  estimating  9  will  allow  continuous  updating  of  the  parame¬ 
ters.  To  achieve  this,  it  is  necessary  to  relate  the  parameters  estimated  at  time  n  to  those 
already  estimated  at  time  n  —  1.  First,  define  the  matrix 


P[n]  =  \4 

i^[n](^[n]]  , 

and  note  from  equation  (A7)  that 

r  <^"’[1]  1 

0[n  -  1]  == 

^T[2] 

5 

ip^[n  —  1] 

so  that 

r  <^'^[1]  1 

P[n  -  1]  =  \ip^{l)ip^[2] . , 

1]] 

<^^[2] 

_  ip^[n  -  1]  _ 

(A12) 


(A13) 


(A14) 
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then 

P“^[n]  =  P~^[n  -  1]  +  ip[n]^'^{n]. 

Recall  from  equation  (All)  that 


6[n] 


y 


(j) 

Y^(p\i]cp^\i] 


U=1 


-1 


=  P[n] 

=  P[n] 


i=l 

n—1 

Y1 


li-l 


v’Wy  w 

Li=l 


+  ip[n]y[n], 


the  third  line  of  the  above  equation  implies 

§[n  —  1]  =  P[n  —  1] 


n—1 


Lt=l 


or 


n—1 

E 

i=l 


(A15) 


(A16) 


(A17) 


(A18) 


=  |^P“^[n]  —  —  !]• 

Where  the  last  line  is  as  a  result  of  using  P-^[n  -  1]  =  P'^[n]  -  ip{n]ip^[n]  from  equa¬ 
tion  (A15).  Using  this  in  the  last  line  of  (A16) 


0[n]  =  P[n]  [P-Hn]^[n  -  1 
=  6[n  -  1]  +  P[n](p[n] 


-  (p[n]tp^[n]9[n  -  1]  4-  (7’Ny[n]j 
y[n]  -  (p^[n]e[n  -  1]]  . 


(A19) 


Note  that  in  the  last  line  ip[n]e[n  -  1]  is  just  a  one-step-ahead  prediction  of  y[n]  based 
on  the  parameters  estimated  at  time  n—1.  The  expression  in  the  brackets  is  then  the 
prediction  error  e[n]  =  yM  ~  y[n].  .The  above  equation  can  then  be  re-written  as 

6[n]  =  9[n  —  1]  +  K[n]e[n],  (A20) 

which  is  appropriately  recursive.  It  remains  to  determine  a  recursion  equation  for  P[n]. 
Recall  from  equation  (A15)  that 

P[n]  =  [P[n  -  1]  4-  (/5[n]¥P^[n]j  .  (A21) 

Using  the  matrix  inversion  lemma  [2] 

[A  +  BCD]~^  -A-^b[C~^ +DA-^B^  ^  DA-^,  (A22) 

where  A,  C  and  {C~^  +  DA~^B]  are  required  to  be  non-singular.  Let 

A  =  P"^[n-1] 

B  =  ^[n]  (^23) 

C  —  I  the  identity  matrix 

D  =  (p^[n], 
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then 

P[n]  =  [F-i[n] +v?[n]<^^[n]] 

=  P[n  —  1]  —  P[n  —  l]<^[n]  +  ip^[n]P[n  —  l]</?[n]j  ip^[n]P[n  —  1]. 

Since  K[n]  =  P[n]<^[n],  we  have 

K{n]  =  P[n  —  l]y>[n]  —  P[n  —  IJv^H  [/  +  ip^[n]P[n  —  l]<(0[n]]  ^  ip^\n]P[n  —  l]</5[n] 

=  P[n  —  l]</?[n]  [/—[/  +  tp^[n]P[n  —  l]93[n]]  ip^[n]P[n  —  l]v5[n]]  (A25) 

=  P[n  -  l]¥p[n]  [/  +  (fi^[n]P[n  -  l]v3[n]]  , 

where  the  last  line  is  as  a  result  of  applying  the  matrix  inversion  lemma  to  the  contents 
of  the  outer  brackets  of  the  line  above.  From  equations  (A24)  and  (A25) 

P[n]  =  P[n  —  1]  —  P[n  —  l]<p[rr]  +  ip^[n]P[n  —  l]<pMj  ip^[r{\P[n  —  1] 

I  +  ip^[n]P[n  -  l](p[n]]  , 


K[n]  =  P[n  —  l](p[n] 
it  follows  that 

P[n]  =  [/  —  A’[n](p[n]]  P[n  —  1]. 

The  full  recursive  procedure  for  estimating  the  parameter  vector  6  is  then, 


(A26) 

(A27) 


1.  Update  K[n] 

K[n]  =  P[n  —  l]^[ra]  +  ip^[n]P[n  —  l]i^[n]j 

for  a  single-input  single-output  system,  this  reduces  to 

K\n^  =  F[n  -  l]ip[n] 

1  -f  ip^[n\P[n  —  l]<p[n] 

9[n\  =  6[n  —  1]  +  K[n]  ^[n]  —  (p^[n]6[n  —  l]j 


2.  Update  6[n] 

3.  Update  P[n] 


P[n]  =  [/  —  iiir[n](p[n]]  P[n  —  1] 
The  procedure  is  initialised  be  setting 


P[0]  = 


a  0  ...  0 

0  a  ...  0 

0  0  ...  a 


=  al 


(A28) 

(A29) 

(A30) 

(A31) 

(A32) 


The  above  algorithm  assumes  knowledge  of  the  input  sequence,  which  is  imavailable  for 
a  stochastic  model.  However,  since  the  prediction  error  for  a  stochastic  model  is  analogous 
to  the  Gaussian  distributed  input  sequence,  it  may  be  used  as  an  estimate  of  the  input 
sequence, 

v[n]  —  e[n]  —  y[n]  —  yF[n\9[n  —  1].  (A33) 

To  initialise,  the  error  sequence  is  assumed  zero  for  n  <  0  ie.,  u[n]  =  0.  Additionally,  the 
initial  parameter  vector  is  also  set  to  zero,  0[O]  =  0. 
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